##replicateOriginalFig1andFig2.R

##This file replicates Figures 1 and Figures 2 from the original paper.
##Confidence intervals are added to the original representations.

library(foreign)
library(Zelig)
library(sandwich)
library(MASS)
library(lme4)
library(VGAM)
library(mlogit)

data <- read.dta(file="clare.dta")
head(data)
dim(data)

############################
## Table 1, Model 2 predicted probabilities (Figure 1 from the paper)
############################
table1Model2Vars <- c("cwinit", "dyad", "rivalry_count1", "weakmr_decline9", "nriv_weakdec9mr", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, table1Model2Vars]
dataNewNoNA <- na.omit(dataNew)

#highrisk <- rep(0, 10)

#12192
#library(Zelig)
#install.packages("sandwich")
#library(sandwich)
table1Model2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA)
coef(table1Model2)


mat <- estfun(table1Model2)
mat <- na.omit(mat)
dim(mat)
N <- nrow(mat)
u <- apply(mat, 2, function(x) tapply(x, dataNewNoNA$dyad, sum))
u <- na.omit(u)
numberOfClusters <- length(unique(dataNewNoNA$dyad))
df <- (numberOfClusters / (numberOfClusters - 1))*((nrow(mat) - 1)/ (nrow(mat) - table1Model2$rank))
vcovCL <- df*sandwich(table1Model2, meat=crossprod(u)/N)
ses <- coeftest(table1Model2, vcovCL)
ses
sesModel2 <- ses




#library(MASS)
#beta.draws <- mvrnorm(10000, mu = as.matrix(table1Model2$coefficients)[,1], Sigma = summary(table1Model2)$cov.unscaled)
beta.draws <- mvrnorm(10000, mu = as.matrix(table1Model2$coefficients)[,1], Sigma = vcovCL)

highrisk3 <- rep(0, 11)
#highrisk3[2:11] <- highrisk
highrisk3[1] <- 1
#values the original authors used in Stata:
highrisk3[5] <- .09
highrisk3[6] <- .5
highrisk3[7] <- .36
highrisk3[8] <- 1  #scalar h_PolRel=1
highrisk3[9] <- 16
highrisk3[10] <- 684
highrisk3[11] <- 42488 #scalar h_PceYrs3=42488 






highrisk3 <- as.matrix(highrisk3)
dim(t(beta.draws))
#test
#p.ests2 <- pnorm(t(highrisk3)%*%t(beta.draws))
#

NUMBERrivals <- seq(0, 6, 2)
yrsWeak <- 1:10
pEsts0 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts2 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts4 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts6 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)

for(NUMBERrivalsIndex in 1:length(NUMBERrivals)){
  for(j in 1:length(yrsWeak)){
    highriskweak <- highrisk3
    highriskweak[2] <- NUMBERrivals[NUMBERrivalsIndex]
    highriskweak[3] <- 0.9^j
    highriskweak[4] <- NUMBERrivals[NUMBERrivalsIndex] * 0.9^j
    #highriskweak[3] <- yrsWeak[j]
    #highriskweak[4] <- NUMBERrivals * yrsWeak[j]
    if (NUMBERrivals[NUMBERrivalsIndex]==0) {
    	  pEsts0[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	} 
    if (NUMBERrivals[NUMBERrivalsIndex]==2) {
    	  pEsts2[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==4) {
    	  pEsts4[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==6) {
    	  pEsts6[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	}       	
  }
}

#pdf(file="Figure1replication.pdf")
plot(yrsWeak-0.1, apply(pEsts0,2,mean), ylim = c(0,0.2), xlim = c(0.8, 10.5), col="red", main="Predicted Probability of Initiation (Original Figure 1)", ylab="Probability of Initiation", xlab="Years Since Challenger Backed Down")
segments(x0 = yrsWeak-0.1, x1 = yrsWeak-0.1,
y0 = apply(pEsts0, 2, quantile, .025),
y1 = apply(pEsts0, 2, quantile, .975), col="red")

points(x = yrsWeak, y = apply(pEsts2,2,mean), col="blue")
segments(x0 = yrsWeak, x1 = yrsWeak,
y0 = apply(pEsts2, 2, quantile, .025),
y1 = apply(pEsts2, 2, quantile, .975), col="blue")

points(x = yrsWeak+0.1, y = apply(pEsts4,2,mean), col="green")
segments(x0 = yrsWeak+0.1, x1 = yrsWeak+0.1,
y0 = apply(pEsts4, 2, quantile, .025),
y1 = apply(pEsts4, 2, quantile, .975), col="green")

points(x = yrsWeak+0.2, y = apply(pEsts6,2,mean))
segments(x0 = yrsWeak+0.2, x1 = yrsWeak+0.2,
y0 = apply(pEsts6, 2, quantile, .025),
y1 = apply(pEsts6, 2, quantile, .975))

legend("topright", c("0 Other Rivals","2 Other Rivals", "4 Other Rivals", "6 Other Rivals"), bty="n",
col=c("red", "blue", "green", "black"), lwd=c(1, 1, 1, 1), lty=c(1, 1, 1, 1))
#dev.off()



############################
## Table 1, Model 2 first differences (Figure 2 from the paper)
############################
##############Plot first differences
highrisk3 <- as.matrix(highrisk3)

#beta.draws is the same as used above
dim(t(beta.draws))
#test
#p.ests2 <- pnorm(t(highrisk3)%*%t(beta.draws))
#

NUMBERrivals <- c(2, 4, 6)
yrsWeak <- 1:10
pEstsFD02 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEstsFD04 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEstsFD06 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)


#NUMBERdecayRates <- c(0.9, 0.3)



for(NUMBERrivalsIndex in 1:length(NUMBERrivals)){
  for(j in 1:length(yrsWeak)){
    #for 0 rivals:
    highriskweak <- highrisk3
    highriskweak[2] <- 0
    highriskweak[3] <- 0.9^j
    #also show strong:  #in this case, all confidence intervals intersect with 0
    #highriskweak[3] <- 0
    highriskweak[4] <- 0 * 0.9^j

    highriskweak2 <- highrisk3
    highriskweak2[2] <- NUMBERrivals[NUMBERrivalsIndex]
    highriskweak2[3] <- 0.9^j
    highriskweak2[4] <- NUMBERrivals[NUMBERrivalsIndex] * 0.9^j


    #if (NUMBERrivals[NUMBERrivalsIndex]==0) {
    	#  pEsts0[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	#} 
    if (NUMBERrivals[NUMBERrivalsIndex]==2) {
    	  pEstsFD02[,j] <- pnorm(t(highriskweak2)%*%t(beta.draws)) - pnorm(t(highriskweak)%*%t(beta.draws))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==4) {
    	  pEstsFD04[,j] <- pnorm(t(highriskweak2)%*%t(beta.draws)) - pnorm(t(highriskweak)%*%t(beta.draws))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==6) {
    	  pEstsFD06[,j] <- pnorm(t(highriskweak2)%*%t(beta.draws)) - pnorm(t(highriskweak)%*%t(beta.draws))
    	}       	
  }
}

#pdf(file="Figure2replication.pdf")
plot(yrsWeak-0.1, apply(pEstsFD02,2,mean), ylim = c(-0.02,0.1), xlim = c(0.8, 10.5), col="red", main="Difference in the Predicted Probability of Initiation (Original Figure 1)", ylab="Difference in the Probability of Initiation", xlab="Years Since Challenger Backed Down")
segments(x0 = yrsWeak-0.1, x1 = yrsWeak-0.1,
y0 = apply(pEstsFD02, 2, quantile, .025),
y1 = apply(pEstsFD02, 2, quantile, .975), col="red")

points(x = yrsWeak, y = apply(pEstsFD04,2,mean), col="blue")
segments(x0 = yrsWeak, x1 = yrsWeak,
y0 = apply(pEstsFD04, 2, quantile, .025),
y1 = apply(pEstsFD04, 2, quantile, .975), col="blue")

points(x = yrsWeak+0.1, y = apply(pEstsFD06,2,mean), col="black")
segments(x0 = yrsWeak+0.1, x1 = yrsWeak+0.1,
y0 = apply(pEstsFD06, 2, quantile, .025),
y1 = apply(pEstsFD06, 2, quantile, .975), col="black")

abline(h = 0, col="orange", lwd=2)

legend("topright", c("2 Other Rivals vs. 0 Other Rivals", "4 Other Rivals vs. 0 Other Rivals", "6 Other Rivals vs. 0 Other Rivals"), bty="n",
col=c("red", "blue", "black"), lwd=c(1, 1, 1), lty=c(1, 1, 1))
#dev.off()